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Abstract. We consider decaying dark matter (DDM) as a resolution to the possible tension 
between cosmic microwave background (CMB) and weak lensing (WL) based determinations 
of the amplitude of matter fluctuations, ug. We perform N-body simulations in a model 
where dark matter decays into dark radiation and develop an accurate fitting formula for 
the non-linear matter power spectrum, which enables us to test the DDM model by the 
combined measurements of CMB, WL and the baryon acoustic oscillation (BAO). We employ 
a Markov chain Monte Carlo analysis to examine the overlap of posterior distributions of the 
cosmological parameters, comparing CMB alone with WL-I-BAO. We find an overlap that 
is significantly larger in the DDM model than in the standard CDM model. This may be 
hinting at DDM, although current data is not constraining enough to unambiguously favour 
a non-zero dark matter decay rate P. From the combined CMBJ-WL data, we obtain a lower 
bound F“^ > 97 Gyr at 95 % C.L, which is less tight than the constraint from CMB alone. 
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1 Introduction 

Observations of temperature and polarization anisotropies in the cosmic microwave back¬ 
ground (CMB) have been a milestone in cosmology, which allow precision tests of various 
cosmological models. The six-parameter A Cold Dark Matter (ACDM) model has emerged 
as the concordance model of the Universe due to its ability to fit the CMB anisotropies and 
various other cosmological observations.^ The recent observation of CMB by Planck allows 
us to determine the cosmological parameters in this model to sub-percent accuracy [1]. 

At the same time, closer analysis of the level of consistency of the CMB and other 
cosmological data has revealed a number of possible tensions in the ACDM model. One 
that has recently received attention is the value of the amplitude of matter fluctuations, 
characterized by the parameter ug. The value of erg estimated from observations of large- 
scale structure such as the weak lensing of distant galaxies and abundance of galaxy clusters 
is significantly smaller than the one inferred from the Planck CMB power spectrum [2-5].^ 
This tension might be indicating a need to move beyond the standard ACDM and to search 
for a mechanism that would suppress the matter fluctuations at late times. 

Several authors have argued that the tension can be alleviated if active neutrinos are 
massive or sterile neutrinos exist [3, 4, 6-8].^ However, this is not the only possible solution for 
the tension, and other models should also be tested on equal footing. A better understanding 
of effects of astrophysical processes might also be important (see, e.g.. Ref. [10]). In this paper 

^ See Ref. [1] and references therein. 

^ Most of these studies were based on the Planck 2013 results, but the discrepancy still exists in the recent 
Planck 2015 data [1]. 

^ The suppression of erg has also been discussed recently in a specific model that combines dark matter 
and dark radiation [9]. 
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we explore another possibility, where dark matter can decay with a lifetime of the order 
of the age of the Universe. 

Such a decaying dark matter (DDM) model is viable if the decay products interact little 
with the visible particles; otherwise the lifetime should be much longer.^ In this paper, we 
restrict ourselves to the simplest class of DDM models, where decay products are invisible 
and massless. In addition, we also assume that all the dark matter consists of DDM. Such a 
model can be characterised with only one extra parameter T in addition to the standard six 
parameters of the ACDM model. While one can consider more complicated models of DDM 
with, e.g., finite masses of decay products or a mixture of DDM and CDM, our analysis of 
the simplest case is sufficient to demonstrate the potential of DDM models. 

In DDM models, the formation of cosmic structure is in general suppressed. This is 
because the decay products free-stream with finite kick-velocity and can escape from the 
gravitational potential wells surrounding matter over-densities. On the other hand, at early 
times before significant decay has occurred, DDM models are completely degenerate with the 
CDM case. This allows the amplitude of matter fluctuations at late time to be suppressed 
relative to the CDM case, without significantly affecting the CMB anisotropy spectrum except 
through the integrated Sachs-Wolfe effect. 

To distinguish DDM from CDM using cosmological observations requires a precise 
knowledge of the matter power spectrum over a wide range of scales. On large scales, struc¬ 
ture formation in DDM models can be described by linear perturbation theory and has been 
studied by many anthors [12-18]. Based on this linear approach and using the Planck CMB 
power spectrum combined with measurements of baryon acoustic oscillations (BAO) and the 
galaxy power spectrum on linear scales, Ref. [18] gives a lower bound T”^ > 160 Gyr on the 
lifetime of DDM.^ 

On the other hand, to make use of other cosmological data requires capturing the 
behaviour of the matter power spectrum to sufficient accuracy on non-linear scales, which 
can be achieved through using N-body simulations.® In this paper we modify the Gadget2 
code [23, 24] to incorporate the effects of DDM. 

N-body simulations of DDM models have recently been used to study the effects on 
halo inner profiles [25, 27], the Lyman-a forest [26] and effects on non-linear matter power 
spectrum [27]. The implications for small-scale problems in the CDM models have also been 
discussed in Ref. [29]. In these stndies, decay products were assumed to be non-relativistic. 
On the other hand, we focus on models with massless decay product, and follows the method¬ 
ology developed in Ref. [28]. We will show that our linear calculations and N-body simula¬ 
tions agree well on scales where non-linearity is subdominant. We develop a fitting formula 
which can describe non-linear matter power spectrum well, using which we are able to obtain 
parameter constraints in the DDM model from observational data, especially including the 
recent CFHTLens weak lensing measurement [30]. 

We note that the possibility for DDM to resolve the tension in the erg estimations has 
also been discussed in Refs. [17, 27, 31], although the models and parameter regions of DDM 
they consider differ from ours. Also, in this paper we focus on the weak lensing measurement, 
which has not been considered in previous DDM studies. 

The structure of this paper is as follows. In Section 2, we define the DDM model and 
discuss the cosmological linear perturbation theory. In Section 3, we describe our N-body 

^For a recent constraint, we refer to Ref. [11]. 

® Ref. [19] gives a comparable constraint, F”^ > 100 Gyr, using only information on distances. 

® For semi-analytic studies based on the halo model see, e.g., Refs [15, 20-22]. 
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simulations and discuss the effects of DDM on the non-linear matter power spectrum. Using 
the simulation results, we obtain constraints on the DDM model from current cosmological 
observations in Section 4. In particular, we examine to what extent the DDM model can 
mitigate the tension in the ug estimations. We conclude in the Section 5. In Appendix A, 
we describe the approximations used in the Boltzmann hierarchy in the linear calculation. 
Appendix B contains details of our fitting formula for the non-linear matter power spectrum. 

Throughout this paper, we denote a derivative with respect to the conformal time r by 
a dot, i.e. '= djdT. We will assume a flat Universe and adiabatic initial perturbations with 
a power-law primordial power spectrum. 

2 Model and linear perturbation calculation 

2.1 Model of decaying dark matter 

We consider a model of dark matter (DM) decaying into dark radiation (DR) with decay 
rate T. We assume the DM is non-relativistic. Note that as long as the decay products are 
massless, the number of DR particles generated from decay of single DM is irrelevant to the 
evolution of the cosmological perturbations. 

Roughly speaking, current observations constrain the lifetime to be larger than 
100 Gyr [18, 19], so we focus on values of around 100 Gyr in obtaining observational 
constraints. 

2.2 Background evolution 

Given a model of DDM, energy conservation for the DM and DR fluids in the unperturbed 
background metric yields 


Pdm T ^T~(-Pdm — 0,^Pdm,i (^-1) 

Pdr T ^T~Lpdr — oJ^Pd m ; (2.2) 


where pi is the unperturbed energy density of fluid i and 'R = d/a is the conformal Hubble 
parameter. 

Due to the energy transfer between DM and DR fluids, pdm and pdr cannot be repre¬ 
sented as explicit functions of a. This means that the present value of the Hubble parameter, 
Hq =100h km/sec/Mpc, cannot be given as an input parameter but is rather a derived pa¬ 
rameter obtained by solving the differential equations (2.1) and (2.2) in conjunction with the 
Friedmann equation 

, (2.3) 

i 

where G is Newton’s constant. 

To specify the background evolution, it is convenient to introduce the following quantity: 


d)j(a) 


Pcrit/h^ 


(2.4) 


where pcrit/h^ = 3(100km/sec/Mpc)^/87rG and Wi is the (constant) equation of state param¬ 
eter of fluid i. Given the initial values uJi{a = 0) = cuj, we can solve the differential equations 
(2.1) and (2.2) for iPdm and ojdr to obtain = Cji{a = 1) and h. Note that if fluid i 
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does not interact with other fluids (i.e., if the r.h.s. in Eqs. (2.1) and (2.2) vanish), Wj(a) is 
constant in time and oji = U.ih?‘. 

The background cosmology is specified by the parameters 

{u)bi ^dmi ^Ai r), (^•5) 

where the subscripts b and A indicate respectively baryons and dark energy. We assume 
that initially LOdr = 0 so that DR is produced only through the decay of DM. Given these 
parameters, Q.dmh'^, ^drh^ and h = can be obtained as derived parameters. For 

later convenience, we define the reduced Hubble constant in the limit of no decay T = 0, to 
be /i 0 = \/Yl,i la this section and in Section 3, we adopt a fiducial set of initial density 
parameters: 

(wft = 0.02216, ujdm = 0.119, WA = 0.318), (2.6) 

and assume an amplitude Ag = 2.43 x 10“® and spectral index n* = 0.963 at a pivot scale 
k = 0.05 Mpc“^ for the primordial power spectrum. This parameter set gives = 0.677 and 
corresponds to the best-fit ACDM model from the Planck CMB temperature power spectrum 
combined with measurements of small-scale CMB power spectrum and BAG [2]. 

In closing this section, in order to avoid confusion, let us clarify the difference between 
uJdm and Qdmh"^- As we noted, the density parameter ojdm{t) is not constant and evolves 
with the time due to the decay of DM. We have defined the initial and final values of uJdm 
as u!dm = ^dm{o- = 0 ) and ^dm.h'^ = w(a = !)• In the limit of no decay T = 0 , codm and 
^dmh'^ coincide. Let us also remind the readers of the definitions of /i 0 = 
h = \/Yli The former is the reduced bubble parameter that are realized in the limit of 

no decay, while the latter is the actual reduced bubble parameter. 


2.3 Linear perturbation evolution and matter power spectrum 

Let us now discuss cosmological perturbation at linear order. We follow the notation of 
Ref. [32] and treat the perturbed quantities in the synchronous gauge of DM. In this gauge, 
the scalar part of the perturbed metric is 

= a(r)^ [dr^ -|- {dij -|- hij(x, r)} dx'^dx^] , (2.7) 

where hij is given in Fourier space as 

hij{k,T) = hi^{k,T)kikj + 6riT{k,T){kikj - (2.8) 

and''denotes a unit vector, i.e. k = k/k. In the following, we focus on a single Fourier mode 
and dependences on k will be abbreviated. 

In this gauge, the continuity equation for DM is 


^dm 


hh 
2 ’ 


(2.9) 


where 6dm = ^Pdm/pdm is the fractional overdensity of DM. Note that Fq. (2.9) is the same 
as in CDM. 

On the other hand, since DR is massless, its perturbation equations can be given in 
terms of its brightness function defined as 


X{h,T) 


1 

Aa^Pdr 


q^dq 

27r^ 


qSfdriq,T), 
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Figure 1. Power spectrum of non-relativistic matter today (z = 0), for DDM models with r[Gyr“^]=3 
(red), 30 (green) and 300 (blue). For reference, the case with the standard CDM is also plotted 
(magenta). 


where 6fdr{Q, a) is the perturbation in the phase-space distribution of DR and h is the unit 
vector of a comoving momentum q = qn. The Boltzmann equation for DR can be given in 
terms of X{h,T) as [14] 

X + ik^X = - ^{iiL + 67)T)At^| (^X - , (2.11) 

where ^ = k ■ n. After multipole expansion, Eq. (2.11) is recast into 

^ ^ (vi - %*o) . (2.12) 

where X(h, r) = with Pi being the Legendre polynomial. 

To solve the perturbation equations, we have modified the Boltzmann code CAME [33]. 
We truncated the Boltzmann hierarchy at £max = 50 up to moderately sub-horizon scales, 
kr < iraa.x- However, there is a difficulty in solving the Boltzmann hierarchy Eq. (2.12) on 
deep sub-horizon scales. A naive truncation of the hierarchy yields fictitious reflections at 
^max, and the accumulation of them quickly blows up Xi. To avoid this, we adopted the 
approximations developed in Ref. [34], details of which are provided in Appendix A. 

In Fig. 1, we plot the matter power spectrum P{k) at present. It is clearly seen that 
as T is increased, P{k) is suppressed on small scales. We have also checked that our code 
reproduces the CMB power spectrum shown in Ref. [14] and the matter power spectrum 
in Ref. [18]. By adjusting the accuracy settings in CAME, we have determined that within 
the parameter range of interest the numerical accuracy of the linear calculation code is sub- 
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percent for the CMB temperature and polarization power spectra, and ~ 1 % for the matter 
power spectrum. 

3 N-body simulation and non-linear matter power spectrum 

3.1 Simulation 

We perform N-body simulations of collisionless particles by using the public code Gadget2 [23, 
24]. To incorporate effects of DDM, we make two modifications, which are also implemented 
previously in Ref. [28]. One is that the mass of N-body particles m is made time-dependent 
and varies as 

m{t)=mi{{l-rdrn) + rdrne~^^}, (3.1) 

where rui is the initial particle mass and rdm = + ^dm) is the initial fraction of the 

non-relativistic matter in DDM. The other is that the expansion of the Universe is given by 
solving Eqs. (2.1) and (2.2). 

In this treatment, effects of perturbations in DR are neglected. However, since DR is 
relativistic and free-streams, on sub-horizon scales — where Newtonian gravity and thus the 
N-body simulation are valid — DR does not cluster and is highly homogeneous and isotropic. 
On the other hand, on super-horizon scales, the perturbations in DR cannot be neglected. 
However, on these scales, perturbations are linear and their evolution can be accurately 
computed within the perturbation theory treatment. Thus, by interpolating results of N- 
body simulations on small scales and linear perturbation theory on large scales, we can 
predict the matter power spectrum in the DDM model over a broad range of scales. 

The initial matter power spectrum at redshift z = 49 is computed using CAME, from 
which initial distributions of N-body particles are generated using the 2LPTic code [35], 
which uses second-order Lagrangian perturbation theory. The simulations contain 1024^ 
particles, and are performed with three different box sizes, L = 1250, 500 and 200 h^^Mpc. 
For each box size, we generate three different realizations, and for each realization of the 
initial conditions, we run simulations with different values of T, which allows us to extract 
effects of the dark matter decay with statistical fluctuations minimised. 

3.2 Non-linear matter power spectrum 

Fig. 2 shows the ratio of the matter power spectrum in the DDM model with the lifetimes 
T”^ = 31.6, 100 and 316 Gyr to that in the CDM model at redshifts z = 0 and 1. The 
power spectrum of N-body simulations is computed using the ComputePK code [38]^ with 
the count-in-cell scheme, and the linear perturbation calculations are performed by using a 
modified version of the CAME code as described in Section 2.3. Note that here we only take 
into account the overdensity of non-relativistic matter (i.e., baryons and DM), so DR is not 
included in computing the fractional overdensity. We also note that the units of horizontal 
axis are Mpc“^, which is different from the ordinary convention of /iMpc“^. This allows us 
to extract the difference in growth of density fluctuation, after removing the effects of the 
background expansion. 

From Fig. 2, we can see that the ratio from the linear calculation 

and that from the N-body simulations coincide in the range k ~ 0.01 — 0.1 Mpc“^, where 
the both the linear and Newtonian approximations should be valid. Since the effects of 

^ http://aramis.obspm.fr/'lhuillier/codes.php 
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Figure 2. Ratio of non-linear matter power spectrum in the DDM model to that in the CDM model. 
Red, green and blue points with error bars are obtained from N-body simulations with box sizes 
L = 1250, 500 and 200 /i^^Mpc, respectively. The dashed magenta line is obtained from the linear 
perturbation calculation in the synchronous gauge by using GAME. The dotted orange line is obtained 
by applying our fitting formula in Appendix B, which exactly overlaps with magenta line on large 
scales. The solid grey horizontal line indicates unity. 


perturbations in the massless decay products become less relevant at smaller scales, this 
agreement validates our treatment of the decay products in the N-body simulations. 
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The ratio p(DDM)/p(CDM) in the linear calculation tends to increase as k decreases 
below 10“^ Mpc“^. This is because the decay product does not free-stream until perturbation 
scales cross the horizon, and its density perturbation more or less traces that of other non- 
relativistic matter (see, e.g., Ref. [17] for detailed discussion). This effect is omitted in the 
N-body simulations, which results in no increase in the ratio on large scales. 

As the scales become smaller, the non-linear matter power spectrum shows a larger 
deviation of the DDM model from the CDM case, while for the linear calculation the ratio 
stays almost constant. This is because the non-linear matter power spectrum is a non-linear 
function of linear power spectrum (with mode-mode couplings) and hence small deviations 
in linear perturbations are in general enhanced in the non-linear power spectrum. As non¬ 
linearity becomes more prominent on smaller scales, the deviation becomes more enhanced. 

It is also apparent that, except for changes in the overall magnitude, the fractional 
deviation p(DDM) j p(CDM) _ almost the same scale-dependence for various values of 

T”^ and z. This allows us to obtain a fitting formula for the non-linear matter power spectrum 
in the DDM model in Appendix B, which is sufficiently accurate for a wide range of T”^ and 
z. This fitting formula is shown by the orange lines in Fig. 2. Our formula can reproduce 
the fractional deviation pfDDM) j p{CDM) _ ^ from the N-body results to within an accuracy 
of 10 %. 

One caveat is that our fitting formula is not tested against cosmological parameters 
other than the fiducial set described above. However, as shown in Fig. 5 in Appendix B, 
the ratio enon-iinear/eiinear (^)5 where e{k) = 1 — p(DDM) ^p(CDM) ^ does not Vary much with 
redshift or F. This suggests that our fitting formula will also depend only weakly on cos¬ 
mological parameters. In addition, for the cosmic shear power spectrum, which measures 
the gravitational potential but not the fractional overdensity of matter, the impact of un¬ 
certainties in the fitting formula would further decrease. This is because deviations of the 
DDM model from the CDM case also arise from the decrease in mean energy density, which 
suppresses the shear power spectra additionally by ~ {(1 — r^^) + ~ 1 — 2rrfmFt. 

The suppression from this effect is in general larger than the contribution from the power 
spectrum of fractional matter overdensity. Thus the weak dependence of our fitting formula 
on cosmological parameters should not significantly affect the following analysis. 

4 Constraints from CMB and weak leasing measnrements 

Having obtained a fitting formula for the non-linear matter power spectrum in the DDM 
model, we now compare the model predictions to the data on weak leasing, CMB and BAO, 
using a modified CosmoMC code [36]. To obtain the non-linear power spectrum in the DDM 
model, we combine our fitting formula for the y'p(CDM) j.g^jQ -^j^h the improved 

HALOFIT formula of Ref. [45] . 

For the CMB data, we use the Planck 2013 temperature maps [39] combined with the 
WMAP 9-year polarization data [40], and denote the combination as CMB. For weak leas¬ 
ing data, we use the tomographic analysis of the CFHTLens data by Ref. [30], which we 
denote as WL. In particular, we adopt the tomographic cosmic shear power spectrum from 
blue galaxy samples with Nt = 6 redshift bins and the “conservative” cut on small angles to 
mitigate systematic errors from baryonic effects. We refer the reader to Ref. [30] for more 
details. We also use the measurements of baryon acoustic scales from galaxy surveys [41-43] 
and the Planck CMB lensing power spectrum [44], which are denoted as BAO and lensing, 
respectively. We use four different combinations of these data, which are WL-I-BAO, CMB 
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alone, CMB+WL, and CMB+WL+BAO+lensing. For the WL+BAO combination, we im¬ 
pose Gaussian priors iVb = 0.0223 ± 0.0009 and Us = 0.96 ± 0.02. Here the ranges for the 
Gaussian priors are taken to be somewhat broader than the lo" errors in the Planck measure¬ 
ments of these quantities. However, the prior range does not affect our final results much 
since WL and BAG are not sensitive to ujb and n^. (In fact, even if we fix the values of ujb 
and Ug in obtaining the constraints from WL-I-BAO, the results do not change much.) 

The cosmic shear power spectrum is obtained as follows. Using the Limber approxima¬ 
tion, the cross-angular power spectrum of the convergence fields k at redshifts zi and Z 2 at 
multipole £ is given as (see, e.g.. Ref. [37]) 

Ci{zi,Z 2 ) = J dr (^l - [47rGa^p{a)]‘^ P{k = ^,z), (4.1) 

where r{z) is the comoving distance to redshift z = l/o — 1, and rj = r(zi), i = 1,2. Note 
that we have omitted the contribution from perturbations of the decay products, because 
current observations are sensitive to perturbations at sub-horizon scales, where those of 
decay products can be neglected. By applying the fitting formula for the non-linear matter 
power spectrum in the DDM model, Eq. (4.1) can provide the theoretical prediction for weak 
leasing observations. 

In order to avoid confusion, let us remind the reader of the definitions of parameters 
which we constrain in the following. oJdm is the initial density parameter of dark matter 
and in general differs from one at present, o's is computed from the matter power 

spectrum of non-relativistic matter; Hq = lOO/i is the actual Hubble parameter at present 
and should not be confused with lOO/i 0 . The definitions of these parameters are presented in 
Sec. 2 . 2 . 

Let us take a look at the parameter constraints in the CDM model first. Fig. 3 shows 
the constraints on parameters of particular interest, ujdm, o'g and ffg- In the codm-o's and 
HQ-ag planes, the posterior distributions from the GMB and WL-I-BAO datasets overlap 
only marginally at the 95 % G.L. level. In particular, the tension resides in estimation of erg. 
This result highlights the existence of the tension between Planck and GFHTLens, which has 
been argued in the literature [1, 2, 5, 7]. To be more quantitative, we explore the posterior 
distributions from these two datasets in the 3-dimensional parameter space of {oJdm^ ctsj 
Hq). We find that the posterior distributions overlap only at more than 90 % G.L., which is 
broadly consistent with the results of Ref. [7], where the full 6 -dimensional parameter space 
is analysed. 

Fig. 4 shows parameter constraints in the DDM model, which has an extra parameter 
F. Compared to the CDM case, the WL-I-BAO dataset allows a larger value of ag. This 
leads to a substantial reduction of the tension between Planck and GFHTLens in the CDM 
model. In the same manner as in the CDM model, we examine the overlap of the posterior 
distributions from CMB and WL-I-BAO in the 4-dimemsional parameter space (F, uJdm, ctsj 
Hq) of the DDM model, and find that they overlap at 55 % C.L. This level of improvement 
appears somewhat better than that achieved by the addition of massive (sterile) neutrinos [7] , 
which reports that posterior distributions from GFHTLens and CMB overlap at 90 % C.L. 
in the massive neutrino model, that has one extra parameter, while they do so at 64 % C.L. 
in the sterile neutrino model, that has two extra parameters. However, a direct comparison 
is complicated by the fact that the analysis uses the shear power spectrum data on smaller 
angular scales than we do, while additionally accounting for theoretical uncertainties due to 
intrinsic alignment of galaxies and the effects of AGN feedback. 
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^dm 


Figure 3. Constraints on the CDM model. The Id and 2d posterior distributions for parameters 
Wdrn, Hq and erg are shown, with the 2d constraints given at 68 % and 95 % C.L. Green, grey, red 
and blue lines correspond to the data sets WL, CMB, CMB+WL and CMB+WL+BAO+lensing, 
respectively. 


We would like to comment on the degeneracy between ag and T from the WL+BAO 
dataset. As seen from the bottom-left panel in Fig. 4, these two parameters are positively 
correlated, which may at first look counterintuitive as matter fluctuations should be more 
suppressed as F becomes larger. However, ug represents the amplitude of the overdensity in 
non-relativistic matter (baryon-l-DM), which is less suppressed than that in the total matter 
(i.e. baryon-|-DM-|-DR) or gravitational potential. Therefore, as F is increased, ag should be 
also increased to fit the amplitude of the cosmic shear power spectrum. 

Finally, we present constraints on the lifetime of DDM, F“^, from the various datasets. 
From the CMB alone, we obtain F~^ > 140 Gyr at 95 % C.L. This is very similar to the 
constraints in Ref. [18], where the authors use CMB-I-BAO combined with the galaxy power 
spectrum. Combining CMB with WL, we obtain a less stringent constraint, F~^ > 97 Gyr. 
This is because the CMB-I-WL data prefers a nonzero value of F, in order to mitigate the 
ag tension in the CDM model. Inclusion of BAG and the CMB lensing data changes the 
lower bound only by a few percent. On the other hand, WL data alone is not yet sufficiently 
constraining even when it is combined with BAO, so we cannot obtain meaningful constraints 
on F“^ above F~^ > 31.6 Gyr, where we can guarantee the accuracy of our fitting formula. 


- 10 - 















0.24 


0.18 


0.12 


78 


it! 


o 72 


66 


0.90 


0.75 


0.60 



0.008 0.016 0.024 0.12 0.18 0.24 66 72 78 0.60 0.75 0.90 


r[Gyr 


Figure 4. The same as in Fig. 3, but for the DDM model, which has an additional parameter F. 


5 Conclusion 

We have studied cosmological structure formation in a model of dark matter decaying into 
massless invisible particles (“dark radiation”), denoted as DDM model. Employing N-body 
simulations, we have obtained an accurate fitting formula for the non-linear matter power 
spectrum in the DDM model. Using this fitting formula, we compared predictions of DDM 
to the tomographic cosmic shear power spectrum measurements from the CFHTLens survey 
in conjunction with CMB and BAO data. We have demonstrated that the tension present 
in estimations of erg in the CDM model is substantially alleviated in the DDM model. 

Quantitatively, we find that the posterior distributions from two different datasets, CMB 
alone and WL-I-BAO, overlap with each other at 55 % C.L. in the DDM model, while they 
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do so only at 90 % C.L. in the CDM model. We thus suggest the discrepancy may be 
hinting at the possibility that dark matter could be unstable. On the other hand, current 
cosmological observations do not provide powerful enough constraints to nail down the finite 
lifetime of DM, Combining the CFHTLens weak lensing data with the Planck and 

WMAP CMB data, we obtain a lower bound, P > 97 Gyr. This is less stringent than the 
one from CMB alone, P > 140 Gyr, which reflects the preference for a non-zero dark matter 
decay width from CMB-I-WL. 

In the near future weak lensing measurements such as the Subaru Hyper Suprime- 
Cam [46], the Dark Energy Survey [47], and Euclid [48], should increase the statistics of 
the data by orders of magnitude. The DDM scenario we have considered in this paper will 
thus definitely be tested by these future observations. One should also point out that certain 
cosmological observations other than weak lensing may also show signs of tension with the 
CMB and the ACDM model. In particular, the number of galaxy clusters observed by the 
Sunyaev-Zeldovich effect and X-ray observations are known to be less than the prediction of 
the best-fit ACDM model to the CMB power spectrum [2]. Since the amplitude of matter 
fluctuations are suppressed in the DDM model, the abundance of dark matter haloes and 
hence hosted galaxy clusters also becomes less than in the CDM model. It would be inter¬ 
esting to consider whether the DDM model can in a consistent manner resolve the possible 
disagreements between all such observations. We leave this discussion to a future work. 
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A Approximations in solving Boltzmann hierarchy 

In this appendix, we summarise the approximations used in solving the Boltzmann hierar¬ 
chy of Eq. (2.12) in the linear perturbation calculation. Our approximations are based on 
Ref. [34]. 

Eirst of all, let us define a new variable X = oJdrX. Then Eq. (2.II) can be rewritten 
as 

X + ik^x = Uldr |?7T -^{hL + 67)t)I^^| + (A.I) 

Noting that Xi vanishes at r 

A(t) 

where 


0, the formal solution of this equation can be given as 
r {A{t') + B(T')^?] , (A.2) 


^drVT “ 1 “ ^dr ^ ; 


B = - 


^dr 


hL + 6 ? 7 t ) ■ 


(A.3) 

(A.4) 
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Since B vanishes at r = 0, after integrating by part and multipole expansion, we obtain 


1 


k 


Xi{t) - B{t)6io + 




■^Jz(A:r) + 


dT'ji{k{T - t')) A{t') - 


B{r') 

A;2 


(A.5) 


where ji is the first kind of spherical Bessel function. 

After horizon crossing kr ~ 1, Xi starts to oscillate as the free-streaming generates 
higher multipole moments from lower ones. In this regime, the ultra-relativistic fluid approx¬ 
imation (UFA) developed in Ref. [34] can be applied. Following Ref. [34], let us introduce a 
new variable Yi = Xi — -^{Bdio + ^BSn), whose behaviour is quite close to that of ji{kT), 
which satisfies j[(x) — ji-i{x) -\- (x) = 0. From Eq. (A.5), one can derive the following 

equation: 


Yi - kYi + —Yi 
r 

= 


/ + 1 


/ dr' — 

'o T-r 


-jl{k{T-T'))\A{T')- 


B{r') 

k^ 


. (A.6) 


The integrand on the r.h.s becomes large only at kr — I < kr' < kr. This allows one to 
approximate the third term on the r.h.s by 


-(/ + !) A(t)- 



0 a: 


(A.7) 


In particular, for I = 2, Eq. (A.6) can be approximately rewritten in terms of Xi as 
^ = kXi - 


^dr 

dVdr J 




6 


dVdr 4 


(A.8) 


where we have assumed \B/k‘^\ \B/{krf\<^\A\. 

On the other hand, deeper inside horizon, perturbations in Xi become very small and 
their effects on the perturbation evolution of non-relativistic matter are not relevant unless the 
energy density of the decay product dominates the Universe. On such scales, the oscillation 
period of Xi is much shorter than the typical time-scale of perturbation evolution of non- 
relativistic matter, so that only non-oscillatory components of Xi suffice in computing the 
matter perturbation evolution. In this regime, the radiation streaming approximation (RSA) 
in Ref. [34] can be applied. Let us consider the Boltzmann hierarchy for Xp 


j ^ ~ hi ^ 5 dm 

^0 — fCAi ^dr ^ \ ^dr . ? 

D 4 

= Ixo, 

where we have assumed Xi = 0 for I >2. Combining these two equations we obtain 

Id'' 


Xn = 


Cbdrk"^ dr 


h-i 4 35dm 
^dr~^ r ^dr ^ 


(A.9) 

(A.IO) 

(A.ll) 
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where we have assumed |Xo| <C A:^|Xo| since we are interested in the non-oscillatory compo¬ 
nents of Xi. By integrating Eq. (A.9), we then obtain 




^dr ^dm 

6k OJdr 4A: ’ 


(A. 12 ) 


where we have used the fact that Xi, ojdrhh and oJdr^dm vanish at r = 0. To compute Hl, we 
adopt the same procedure as in Ref. [34], which is valid in a matter-dominated Universe. 

In our calculation, the Boltzmann hierarchy Eq. (2.12) truncated at I = /max is directly 
solved from superhorizon scales up to kr = /max- Then we switch the UFA on and the 
Boltzmann hierarchy is truncated at / = 2 by adopting Eq. (A. 8 ). Finally, when kr > 10 x/max 
in the matter dominated Universe, we switch the RSA on and substitute Eqs. (A.11)-(A.12) in 
the source term of the Einstein equation. These approximations avoid the unwanted blowup 
of Xi caused by the truncation of Boltzmann hierarchy. Note that to employ the RSA, we 
need to assume that the decay product is a minor component of the energy density of the 
Universe. This assumption is valid in the parameter region we explore in this paper. 


B Fitting formula for non-linear matter power spectrum in DDM model 

In this appendix we present the fitting formula for the non-linear power spectrum of non- 
relativistic matter in the DDM model. First, let us present the fitting formula for the 
suppression in the linear power spectrum Riinpar at the small scale limit. We find that 
eiinear = 1 " ^iSnear^V^u^f ^ at fe ^ oo Can be approximately given as 

eiinear(r,z) —) , (B.l) 

where a, /3 and 7 are functions of Ub, /i 0 and Um = + ^dm in the case of our model, which 

can be fitted by 

a{ujb, = (5.323 - 1.4644u - 1.391u) -h (-2.055 -h 1.329rt 0.8672u)u; 

-h(0.2682 -0.3509ri)t(;2, (B.2) 

I3{ujb, h,ujm) = 0.9260 (0.05735 - 0.02690u)r(; (-0.01373 0.006713u)t(;^ (B.3) 

'y{u}b,h,ujm) = (1.653- 0.7860u)-h (0.4884-F0.1754 u)u) 

-^(-0.2512-^0.07558^;)^c^ (B.4) 

where u = cj^/0.02216, v = /i0/O.6777 and w = a;m/0.1412. The accuracy is 10 % for 
T”^ > 31Gyr and z < 1 as long as 0.019 < ujb ^ 0.026, 0.6 < /i 0 < 0.8 and 0.09 < uim < 0.28. 

In Fig. 5 we show the suppression in the non-linear power spectrum enon-iinear(/>:) = 
1 - ^n™Lr/^norI?MLr “ units of cunear for T"! = 31, 100, 316 Gyr and z =0, 1. The 
cosmological parameters in Eq. (2.6) are adopted here. One can see that enon-iinear(/i:)/eiinear 
has only a weak dependence both on T”^ and 2 ;. We find enon-iinear(fc)/eiinear can be fitted by 
the following functional form: 


enon-linear(/s) _ l + a{k/MpC 


^linear 


1 -|- 6 (A:/Mpc ^)'? ’ 


(B.5) 








Figure 5. Enhancement of suppression from linear to non-linear matter power spectrum 
enon-iinear(fc)/eiinear for =31.6, 100, 316 Gyr and z =0, 1. Here we adopt the result of N-body 
simulations with a box size L = 200 h^^Mpc. 


where a, 6, p and q are functions of F and z as 


a(r, z) = 0.7208 + 2,027 (^) + 3-«l ( 

b{T, z) = 0.0120 -h 2.786 ( ) -h 0.6499 ( —*— 

VGyr V \l + z 

p{T, z) = 1.045 + 1.225 + 0.2207 

g(r, z) = 0.9922 + 1.735 + 0.2154 


(B.6) 

(B.7) 

(B.8) 

(B.9) 


The accuracy of this fitting formula is better than 10 % for F”^ > 31 Gyr and z < 1. We 
note that this fitting formula is calibrated with simulations only for the fiducial cosmological 
parameters in Eq. (2.6). However, having observed the fact that enon-iinear(^)/eiinear varies 
only weakly for various F~^ and z, we speculate that this quantity does not vary much for 
different cosmological parameters, as long as the linear power spectrum does not deviate 
significantly from the fiducial one. This condition is effectively satisfied when one compares 
model predictions with observed data of matter power spectrum or CMB. 
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